clear
rng default
n_types=4;
all_years=1940:1:1967;
n_years=size(all_years,2);

%% Construct pictures reporting the bidimensional projections of the identified sets computed in the cluster. 
%When the hitting the grid bounds report -Inf, +Inf in the table below. 
%3 indicates -Inf, +Inf
%2 indicates lower bound -Inf
%1 indicates upper bound +Inf
%NaN indicates finite lower and upper bound

%U	2	3	4	2	3	4	2	3	4	2	3	4
U=[ 1	3	3	1	3	3	1	3	3	1	1	1; ...
 	1	3	3	1	3	3	1	3	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	1];

%V	2	3	4	2	3	4	2	3	4	2	3	4
V=[ NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1];


U_reshape=cell(n_years,1);
for h=1:n_years
    U_reshape{h}=NaN(n_types, n_types);
    for i=1:n_types
        U_reshape{h}(i,2:n_types)=U(h,(i-1)*(n_types-1)+1:i*(n_types-1));
    end
end

V_reshape=cell(n_years,1);
for h=1:n_years
    V_reshape{h}=NaN(n_types, n_types);
    for i=1:n_types
        V_reshape{h}(2:n_types,i)=(V(h,(i-1)*(n_types-1)+1:i*(n_types-1))).';
    end
end

%%
fp='.../SpecB/output'; %specify path

C_21_m=zeros(n_years,2); %[min max] for each year
C_32_m=zeros(n_years,2); %[min max] for each year
C_43_m=zeros(n_years,2); %[min max] for each year
C_21_w=zeros(n_years,2); %[min max] for each year
C_32_w=zeros(n_years,2); %[min max] for each year
C_43_w=zeros(n_years,2); %[min max] for each year

for h=1:n_years
    
    IdSetM=cell(n_types,1);
    IdSetW=cell(n_types,1);
    C_interval_m=cell(n_types,n_types);
    C_interval_w=cell(n_types,n_types);
    
    load(num2str(all_years(h)), 'PYX_cond', 'PXY_cond');  
    
    for j=1:n_types
        
       filepath = fullfile( fp, ['IdSetM' num2str(all_years(h)) '_final_' num2str(j) '.mat'] );
       if exist( filepath, 'file' )
          A = getfield( load( filepath ), 'IdSetM_final' ); 
          IdSetM{j}=A;
       end 
       
       filepath = fullfile( fp, ['IdSetW' num2str(all_years(h)) '_final_' num2str(j) '.mat'] );
       if exist( filepath, 'file' )
          A = getfield( load( filepath ), 'IdSetW_final' ); 
          IdSetW{j}=A;
       end  
    end
    
    
    for i=1:n_types-1
    for k=i+1:n_types
        %Men
        U_i=IdSetM{i}(:,2:end); %U_1 U_2 U_3 U_4
        U_k=IdSetM{k}(:,2:end);

        C_tempk=sum(PYX_cond(k,1:end-1).*U_k,2);
        C_tempi=sum(PYX_cond(i,1:end-1).*U_i,2);
        C_interval_m{k,i}=[min(C_tempk)-max(C_tempi) max(C_tempk)-min(C_tempi)];
        
        %Replace +Inf, -Inf
        if any(U_reshape{h}(k,:)==1) || any(U_reshape{h}(k,:)==3)  || any(U_reshape{h}(i,:)==2) || any(U_reshape{h}(i,:)==3) 
           C_interval_m{k,i}(2)=Inf;
        end
       if  any(U_reshape{h}(k,:)==2) || any(U_reshape{h}(k,:)==3)  || any(U_reshape{h}(i,:)==1) || any(U_reshape{h}(i,:)==3) 
           C_interval_m{k,i}(1)=-Inf;
       end
        
        %Women
        V_i=IdSetW{i}(:,2:end);
        V_k=IdSetW{k}(:,2:end);

        C_tempk=sum(PXY_cond(k,1:end-1).*V_k,2);
        C_tempi=sum(PXY_cond(i,1:end-1).*V_i,2);
        C_interval_w{k,i}=[min(C_tempk)-max(C_tempi) max(C_tempk)-min(C_tempi)];
        
        
        %Replace +Inf, -Inf
        if any(V_reshape{h}(:,k)==1) || any(V_reshape{h}(:,k)==3)  || any(V_reshape{h}(:,i)==2) || any(V_reshape{h}(:,i)==3) 
           C_interval_w{k,i}(2)=Inf;
        end
       if  any(V_reshape{h}(:,k)==2) || any(V_reshape{h}(:,k)==3)  || any(V_reshape{h}(:,i)==1) || any(V_reshape{h}(:,i)==3) 
           C_interval_w{k,i}(1)=-Inf;
       end
       
    end
    end
    
    C_21_m(h,1)=C_interval_m{2,1}(1);
    C_21_m(h,2)=C_interval_m{2,1}(2);
    C_32_m(h,1)=C_interval_m{3,2}(1);
    C_32_m(h,2)=C_interval_m{3,2}(2);
    C_43_m(h,1)=C_interval_m{4,3}(1);
    C_43_m(h,2)=C_interval_m{4,3}(2);
    
    C_21_w(h,1)=C_interval_w{2,1}(1);
    C_21_w(h,2)=C_interval_w{2,1}(2);
    C_32_w(h,1)=C_interval_w{3,2}(1);
    C_32_w(h,2)=C_interval_w{3,2}(2);
    C_43_w(h,1)=C_interval_w{4,3}(1);
    C_43_w(h,2)=C_interval_w{4,3}(2);

   
   
end

save('Cm_R2', 'C_21_m', 'C_32_m', 'C_43_m')
save('Cw_R2', 'C_21_w', 'C_32_w', 'C_43_w')

